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Vh ' ABSTRACT 

Oh. 

Aims. The detection of significant concentrations of crystalline silicates in comets indicates an extensive radial mixing in the pri- 
mordial solar nebula, i.e. the protoplanetary accretion disk of our solar system. In studying the radial transport of matter within 
protoplanetary disks by numerical model calculations it is essential to resolve the vertical disk structure since matter is mixed radi- 
ally inward and outward by a complex 2-dimensional flow pattern within the disk that is superposed on the global inward directed 
accretion flow. It is further essential to follow numerically the advection-diffusion processes over a period of at least 10 6 yrs to allow 
' pH ' ' for a full development of the radial concentration profile built up by radial mixing from the warm inner to the cool outer parts of 

protoplanetary disks beyond of 10 AU. 

Methods. Numerical model calculations of protoplanetary accretion disks with radial and vertical mixing are performed by solving a 
set of 2-dimensional transport-diffusion-reaction equations for some important tracers self-consistently with the set of disk equations 
in the 1 + 1-dimensional approximation. The global 2D velocity field of the disk is calculated from an approximate analytical solution 
for the meridional flow pattern, which exhibits an inward drift in the upper layers and an outward drift in the midplane in most parts 
of the disk. The disk model is based on the /^-prescription of viscosity and considers vertical self-gravitation of the disk. This kind of 
semianalytical approximations allows with presently available computer capacity to follow the evolution of the disk and the transport 
and mixing of tracers in vertically resolved 2D-models over the required long periods of disk evolution. The mixing processes in the 
disk are studied for the following species: amorphous silicate grains (forsterite, enstatite) which crystallise by annealing in the warm 
inner parts of the disk, and carbonaceous grains which are destroyed by surface reactions with OH molecules at elevated temperatures. 
Results. Considerable fractions of crystallised silicates and methane (formed as a by-product of carbon combustion) are transported 
to the site of comet formation far from the protosun within a period of 10 6 yrs. The 2-dimensional transport of tracers in the solar neb- 
ula therefore offers a natural explanation for the presence of crystalline silicates in comets and the significant portions of crystalline 
silicates observed in accretion disks around young stellar objects. 
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1 . Introduction matter in protoplanetary disks by advection occurs on a similar 

. £h ' timescale as the transport by turbulent diffusion, namely on the 

k> In the standard one-zone model of protoplanetary accretion disks viscous timeS cale. Hence, besides a realistic description of the 

r £ (Pnngle 1981 ; Lin & Papaloizou |T985) the vertically averaged turbulen t diffusion, an exact knowledge of the large-scale flow 

C3 ; flow field of the disk is characterized by a radial inward drift fieW in disks is essential for ca i cu i ating the transport of matter 

of the disk matter within the inner, chemically active zone of within disks Here we consider the effect of meridional flows. 

the disk. Only in the icy region far away from the proto-sun Another type of hydrodynamic mixing associated with gravita- 

the disk's flow field is directed outward. However, Urpin (1981 tional instabilities is discussed by Boss d2004l 120071 120081 and 

found in his analytical work large-scale meridional flow pat- f ounc [ to be very efficient in that case 

terns to exist. This meridional flow field is characterized by an 

outward directed drift of the disk matter close to the disk mid- p rev ious model calculations of protoplanetary disks which 
plane whereas the flow is directed inward in higher layers. The indude the calcu i ation of radia i mixing of spe cies have consid- 
result of Urpin (1984) has been confirmed by several authors ered only the vert ically averaged one-zone velocity field as the 
by different analytical, semi-analyticaland numerical methods flow field of the disk (St evenson & Lunine 1988 1 Cyr et al. 1998; 
(Siemiginowska 1988 ; Kley & Lin 11992) Rozyczka et al. 11994) Drouart [T999t Bockelee-Morvan et al. 2002 as well as the series 
Kluzmak & KitaUXXS Regev & Gitelman|2002l Tscharnuter & of papers of the ITA group: Gail 2 001 ; Wehrstedt & GaiH2002l 
Gail 120071). Thus a meridional flow field seems to be the urn- Gail 2QQ2 . Ga0 2004; Wehrste dt & Gail 2003 , henceforth called 
versal type of flow pattern existing in protoplanetary accretion p apers T _ V ). The only exception are the 2-dimensional model 
dlsks - calculation of Keller & Gail (120041 henceforth called Paper VI) 
With respect to radial mixing processes in protoplanetary and Tscharnuter & Gail (120071) where for the first time the merid- 
disks the structure of the flow field is of utmost importance. This i ona i velocity field is used for the computation of the radial mix- 
is a consequence of the fact that the average radial transport of i ng f spe cies in protoplanetary disks. The results show that the 
outward transport of species is much more efficient in the merid- 
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In the present work we extend the work of Paper VI. This 
is done by improving the time-dependent one-zone models of 
Papers II and V by calculating the disk structure in the 1 + 1- 
dimensional approximation (e.g. Lin & Papaloizou 1985 ) simul- 
taneously with the 2-dimensional mixing of species in the disk. 
We therefore accept the small deviations of the 1 + 1-dimensional 
disk structure from the disk structure of the exact 2-dimensional 
hydrodynamic calculations of Paper VI in order to save comput- 
ing time. In contrast, however, the present 1 + 1-dimensional disk 
model is coupled with a sophisticated chemical model which in- 
cludes the calculation of equilibrium condensation of the most 
abundant solids, annealing of silicates, and combustion of solid 
carbon. Furthermore, the present disk model includes a detailed 
opacity calculation considering the Rosseland and Planck opac- 
ity means of the most important absorbers in the disk. 

With this powerful tool we investigate radial mixing pro- 
cesses in protoplanetary disks. In particular the radial mixing 
of species in the solar nebula is of great interest to explain the 
composition of the most pristine solar system bodies, i.e., the 
comets. From observations it is long known that some fraction of 
the dust in comets is crystalline (e.g. Swamy et al. 1988; Hanner 
et al. 119971 ) and crystalline silicate grains have been detected 
in accretion disks around young stellar objects (e.g. Meeus et al. 
HJOl] [Bouwman et al. l20(fl1 van Boekel et al. l2TJ04ll2005l Keller 
et al. 2005). Both these observations are interpreted as resulting 
from mixing material from the central parts of the disk to the 
outer regions though also other processes have been invoked for 
that (see, e.g., Alexander et al. 120071 Wooden et al. 120071 for a 
discussion). 

The paper is organized as follows: Section [2] presents the 
treatment of the 2-dimensional transport of tracers in the disk 
within the disk model. In Sects. 0and[4]the method of calculat- 
ing the radial and vertical disk structure is described. Section [5] 
addresses the calculation of the velocity field of the disk, in par- 
ticular the meridional flow field. Section[6]deals with the numer- 
ical treatment of the disk model. In Sect. [7] we present the results 
and finally make our conclusions in Sect. [8] 

2. Transport of tracers 

2.1. The transport-reaction equation 

The transport of tracers with small concentration embedded in 
a carrier medium is given by the transport-reaction equation 
(Hirschfelder et al. [T9P1 cf . Paper VI) 

— riCi + V • Hep, - V ■ nDj Vc,- + R, , (1) 
dt 

where n denotes the particle density of the carrier medium, v, 
the velocity vector of tracer component i, D, the binary diffusion 
coefficient of tracer i relative to the carrier, /?, the rate term of 
gains and losses by chemical reactions of tracer i, and 

et = — (2) 

n 

the concentration of tracer i. rtj denotes the particle density of 
tracer i. The second term on the l.h.s. of Eq. (fl]i is the advection 
term. The first term on the r.h.s. the diffusion term. 
By using the continuity equation 

n 

— + V • nv = (3) 

dt 

the transport-reaction equation (0 can be changed to 

_£i + v . V Ci = _ V ■ nDi Vc,- + — . (4) 

at n n 



For simplicity 



(5) 



is assumed, i.e., all tracer velocities v,- equal the velocity v of 
the carrier medium. For the present model calculations this is an 
acceptable assumption since the considered tracers are micron- 
sized dust grains which are carried along with the flow in most 
parts of the disk. More precisely, for micron-sized particles 
Eq. 01 only holds in the densest and chemically active disk re- 
gions whereas in the less dense outskirts of the disk the motions 
of carrier gas and tracers may decouple, i.e., v, + v. 

The binary diffusion coefficient is assumed to be given by 



D t = — , 

St' 



(6) 



where v is the kinematic viscosity of the disk matter and S , the 
Schmidt number of tracer i. The Schmidt number is set in the 
present calculations to 



S t = l 



(7) 



for all tracers. The approximation for the Schmidt number is 
based on the same approximations as assumption 01. A Schmidt 
number equal to unity means that the tracer moves as the carrier 
gas does. In less dense parts of the disk this approximation fails, 
and S t > 1 . 

Since regions of low density do not contribute much to the 
total tracer transport, we apply (0 and (0 for all tracers in 
the present model. The value of S , can, however, not be fixed 
with any precision because of the unclear physics of the viscous 
transport in accretion disks. For some discussions on the value 
of Si see, e.g., Johansen et a l.|2005l [20061 Turner et al. |2006l 
Pavlyuchenkov & Dullemond 2007 

For calculating the tracer transport in two dimensions we as- 
sume axial symmetry and use polar coordinates (/ , 9). Here, r' 
denotes the polar radius and 9 the polar angle measured from 
the midplane of the disk. With this choice of coordinates the 
transport-reaction equation © takes the form 



dcij ; dcij 1 dctj 1 d , 2 dcij 

■ + v r - + vg— — — —^— t— r nD -^— 



8t 



8r' 

1 



V 89 



2 n dr' 



Br' 



d do/ i R/ i 

— — cos 6 nD — — + — 

r' z n cos 9 89 o9 n 



(8) 



where v' r and vg are the velocities in radial and polar direction, 
respectively. The index j denotes a kind of tracer occurring in a 
special modification i (see following section). 

2.2. The set of tracers 

The tracers which are considered in the model calculations are: 

1 . silicate dust grains (forsterite and enstatite) of different de- 
grees of crystallisation and 

2. carbon dust grains of different sizes. 

The tracers experience the following reactions: 

1 . Silicate dust grains (forsterite and enstatite) start to anneal at 
temperatures above ~ 800 K, i.e., the degree of crystallisa- 
tion increases until it reaches unity. 

2. Carbon dust grains become decomposed by reactions with 
OH molecules at the grain's surfaces. This process critically 
depends on the density of OH molecules in the gas phase and 
starts to operate at temperatures of ~ 1 100K under condi- 
tions encountered in protoplanetary disks. 
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The rate terms for the basic reactions are given in Paper II. 

The set of transport-reaction equations ([8]l determines the 
concentration of crystalline forsterite, Cj,for, and crystalline en- 
statite, c, ens , with different degrees of crystallisation x,, and the 
concentration c, car of solid carbon grains with different grain 
sizes a,. The cy are used to calculate the average degrees of 
crystallisation of forsterite and enstatite (cf. Paper II) 



/cry,for — / t %i Ci,forj ? . Q,for 



i=l 



'Sli 'Sll 

/cry.ens — / i %i Q.ens/ / Q.ens ? 



and the degree of condensation of C into solid carbon 



/car — 



1 



lax- 1 



Vo. car ec 



/=! 



V _ 3 

/ j n a i C i,car 



(9) 



(10) 



(11) 



Vo.car denotes the volume of a carbon atom within solid car- 
bon and ec the (solar) abundance of C. 7 s ji is the number of 
sampling points for the discretised degrees of crystallisation x,-, 
i - 1 . . . / s a, and 7 car the number of sampling points for the size 
spectrum of solid carbon grains a,, i - 1 . . . / car . The set of sam- 
pling points for x, and a, is chosen as in Paper V (I s n = 3, 

/car = 31). 

For a more detailed discussion of the processes of silicate 
annealing and carbon combustion the reader is referred to the 
other papers of this series, in particular Paper I and II. 

2.3. Solution of the set of transport-reaction equations 

The set of transport-reaction equations (O is discretised in first 
order in time and in second order in space with respect to the 
diffusion terms. The advection terms are treated by a standard 
upwind method (see Paper II). 

The transport-reaction equations ([8]) are solved by an ADI 
(Alternating Direction Implicit) method (Press et al. 1992J. By 
ADI each transport-reaction equation is split into two equations 
for separate directions in space: one equation contains the terms 
in r' -direction and the other equation the terms in f?-direction. 
The rate term is split into two half-steps and equally distributed 
on both directions. The resulting matrices of the separate equa- 
tions have tri-diagonal structure and thus can be inverted nu- 
merically fast and easily (Press et al. 1992). Inverting the ma- 
trix of the unsplitted transport-reaction equation ([8j, which has a 
band-tri-diagonal shape, would be numerically much more time- 
consuming and therefore less efficient than the ADI method. 

To ensure numerical stability the sequence of solution of 
the splitted equations is permuted between two successive time 
steps (therefore Alternating Directions Implicit method). The 
transport-reaction equations are solved fully implicit. 

3. Radial disk structure 

The radial structure of the disk is calculated in the one-zone ap- 
proximation by using the /^-prescription of viscosity and by ac- 
counting for the vertical self-gravity of the disk. The resulting 
set of equations for the radial disk structure is given by (cf. Paper 
V): 
1.) Time evolution of the surface density X: 

51 3d d 

-j- = -—ir—vZ^r. (12) 

at r or or 



where r is the radial distance from the mass-centre. 
2.) Keplerian angular velocity: 



Q= — 
r 



GM t 



(13) 



where vk is the Keplerian velocity, G the gravitational constant 
and M, the stellar mass. 
3. ) j3- viscosity: 



Vf) = j3r l Q.. 



(14) 



where f3 denotes the viscosity parameter which is chosen to be 
10~ 5 in the model calculations (cf. Paper V). 
4.) Isothermal sound speed: 



k B T c 



(15) 



where k B is the Boltzmann constant, T c the temperature in the 

midplane of the disk, /u the mean molecular weight and mn the 

proton mass, respectively. 

5.) Pressure scale height by accounting for the disk's vertical 

self-gravitation: 



2nGL 



11 + 



c 8 Q V 

InGL) 



1 



6.) Mean vertically averaged mass density: 

S 
Pm " 2h s • 
7.) Mean molecular weight: 

PmksTc 



M = 



mn(pH + PH 2 + PHe) 



(16) 



(17) 



(18) 



The calculation of the partial pressure px of species X in chem- 
ical equilibrium is described in Paper I. 

8.) Rosseland and Planck means of the mass extinction coeffi- 
cient: 



KR/P,for - /cry.for^R/P.for.cry 

+ (1 _ /cry.forj ^R/P,sil,am 
^R/P,ens — /cry,ens^R/P,ens,cry 

+ ( t — /cry.ens 1 ^R/P,sil,am 
KR/P,dust = /car"'R/P,car + /ens^R/P.ens + /for^R/P.for 
+ /iro'<'R/P,iro + /cor*'R/P,cor 

KR/P - /ice^R/P.ice + (1 -/ice)*"R/P,dust 
' \t /cor/ ^mol ■ 



(19) 

(20) 
(21) 

(22) 



For the calculation of the opacity the main dust absorbers are 
considered in the disk model. The abbreviations used here de- 
note: 'for' = (crystalline) forsterite, 'ens' = (crystalline) en- 
statite, 'sil,am' = amorphous silicate, 'car' = solid carbon, 'iro' 
= solid iron, 'cor' = corundum and 'ice' = water ice, 'mol' = 
molecules, fz denotes the degree of condensation of the key el- 
ement of condensate Z (Si for the silicates, C for solid carbon, 
Fe for solid iron, Al for corundum and O for water ice). The 
method of calculating / car , /; 10 and /i ce is described in Paper I 
and II, respectively, and of ff or , / ens and / cor in Paper V. / cr y,for 
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and /cry.ens are calculated by Eqs. © and ([TOb , respectively, and 
/ car by Eq. ( fTTT i. Analytical fit formulae for the Rosseland and 
Planck mean of the opacity of the individual species are given in 
Wehrstedt ([20031 . 

9.) Rosseland and Planck mean vertical optical depth at the mid- 
plane: 



_ l v 
tr/p - t-^kr/p • 

10.) Viscous dissipation rate: 

E v = - Q 2 vS . 
8 

11.) Effective temperature of the disk surface: 



(tT, 



eft 



1 + 



1 

4n 



E v + crT. 



cloud 



where a is the Stefan-Boltzmann constant and T, 
ature of the ambient molecular cloud. 
12.) Temperature at the midplane: 



(23) 
(24) 

(25) 
,i oud the temper- 



ed = tTr + 



1 



e v + o-n 



[4 4Tpy 

13.) Radial drift velocity: 

3 d v r 
— vLyr, 



14.) Tracer transport in the one-zone approximation: 



~ +Vr ~dr~ 



1 d n daj Rtj 

— — mD — — + — - 
rn or or n 



15.) Rate of mass accretion: 
M(r,t) = 2nrLv r . 



(26) 



(27) 



(28) 



(29) 



16.) Variation of the stellar mass by accretion of matter onto the 
protostar: 



M,(f) = M„(0) - \ M(r in , t') dt' . 
Jo 



(30) 



Note that M(r\ n ) is negative. 

The set of equations for the radial disk structure ([T2l - (f3Qb 
is solved with standard methods up to an accuracy of A ra( j = 10~ 5 
for the radial model structure. 



4. Vertical disk structure 

Since we wish to investigate the radial and vertical transport of 
tracers in protoplanetary disks as well as its feedback on the 
disk structure, it is necessary to resolve the vertical structure of 
the disk. For this purpose we choose the 1 + 1 -dimensional ap- 
proximation for the calculation of the disk structure (e.g. Lin 
& Papaloizou 1985). In this approximation the vertical strati- 
fication is calculated separately for each grid point of the radial 
one-zone model (see Sect.0 by assuming a hydrostatic structure 
in the vertical direction. Therefore we accept some deviation of 
the results of the 1 + 1-dimensional model calculation from that 
of the exact 2-dimensional hydrodynamic calculation, but avoid 
the numerical complexity of the latter. 



4. 1 . Vertical disk equations 

In the following the set of equations for the vertical disk structure 
in the 1 + 1-dimensional approximation at a certain radius r is 
given. 

We define a z-dependent column density: 



o-(z) 



f 



p(z')dz' 



(31) 



where p(z) is the density at the height z above the midplane at 
the given radius r. Note that cr{z) - j^(r) for z — > with this 
definition. The value of E(r) is taken from the one-zone model. 
The differential form of Eq. d3TT > is 



do- 

Tz=~ p - 

The pressure stratification is determined by 



~ = -Q. 2 z-4nG(h-o-(z) 
p oz \2 



(32) 



(33) 



where P denotes the pressure. The second term on the r.h.s. of 
Eq. ( 1331 is the vertical acceleration of a self-gravitating infi- 
nite slab (Paczyriski 1978). The (vertical) self-gravity of the disk 
is taken into account since in Paper V it has been found that 
the self-gravity significantly modifies the disk structure for disk 
masses above ~ 0.01 M,. Since in the present model calculations 
we choose an initial disk mass of Md; s k = 0.2 Af», the self-gravity 
of the disk has to be considered. 

The generation of heat in the disk is assumed to be caused 
exclusively by viscous friction. Thus the vertical gradient of the 
energy flux F z is given by (e.g. Lin & Papaloizou l 19851 1 



dF z 9 2 

1F = 4 QV " 



(34) 



Note that E v - F z (z — * oo) holds according to Eq. d24l . 
Further note that irradiation of the disk by the star is neglected 
in Eq. ([34]). 

The solution of Eq. ( 1341 requires the specification of the ver- 
tical dependency of the viscosity v. However, the mechanism of 
turbulence in disks as the main source of viscosity is a matter of 
much debate so far (see e.g. Richard & Davis 2004; Johansen et 
al. 2005, 2006; Turner et al. 2006; Pavlyuchenkov & Dullemond 
2007 and references therein). Due to the lack of a precise knowl- 
edge of the vertical variation of v, the viscosity at point (r,z) is 
set equal to the value which follows from the one-zone approxi- 
mation at r (Eq. (fT4l). i.e., the viscosity is set vertically constant. 

The heat is assumed to be transported solely by radiation. 
Heat transport by convection is neglected since it is inefficient at 
low temperatures. In this case, the vertical gradient of the tem- 
perature T(z) is given by (e.g. Lin & Papaloizou 1985) 



dT(z) 
dz 



3_ K R pF z 

16 o-T^z) 



(35) 



The Rosseland mean opacity at each vertical location is calcu- 
lated as in the radial model (Eqs. ( fT9b - ([221 ). 

The set of vertical disk equations is closed by the equation 
of state for an ideal gas 



P 



(36) 
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4.2. Solution of the set of vertical disk equations 

The set of equations for the vertical disk structure ( f32b - d36l > 
constitutes a two-point boundary value problem. We define the 
top of the atmosphere at the height zo where the vertical opti- 
cal depth is close to zero. According to the Eddington-Barbier 
approximation (e.g. Mihalas 1978) the temperature at r = is 
determined by 



rjiA 'T'4 . r T<A 

1 ~ 2 eff cloud 



(37) 



The column density at zo has a very small value. In the model 
calculations we assume for <x a fixed value of 



C() 



-■ 1(T 6 I. 



The energy flux at zo follows from Eq. ( l34t to be 



F = -Q 2 v l-S-o-o 



(38) 



(39) 



Finally the pressure at z has to be determined. By assuming an 
isothermal stratification for small optical depths Eq. ( 1331 ) can be 
integrated. Its solution is 



Pn = 



with 



2CqO- 



1 



71 h 1 - <t>(jc ) 



(40) 



2nGcr<jho 1 z , c 

x ° = 2 + "FT" ' h ° = 7=T ' 

ci. V2 h o Q 



co = 



k B T 



<t> denotes the error function and po the mean molecular weight 
at zo which is set to 7/30] 

Starting with an initial guess of the height zo, the set of verti- 
cal disk equations (T32t - d36i l is solved by a Runge-Kutta method 
of 5th order (Press et al. 1992) from zq down to the midplane of 
the disk at z = 0. The Runge-Kutta solver is equipped with a 
step size control which keeps the number of vertical grid points 
small. In general, the column density <x(z) at z = found by 
this procedure is not equal to the midplane value E/2. Then zo 
is iterated by a shooting method until cr(z — 0) = E/2 within an 
accuracy of the vertical model of A vert = 5 • 10" 5 . 

5. Velocity Field 

Finally, the velocity field in the 2-dimensional disk model has to 
be specified. We will show in this paper that the flow field in the 
disk affects the radial mixing of tracers in the disk, exceeding 
the contribution of turbulent diffusion. Therefore the actual flow 
field in the protoplanetary disk has to be determined. 

5.1. Velocity field of the one-zone approximation 

The velocity field of the one-zone approximation is given by 
Eq. $T7\ . The main feature of this velocity field is that there ex- 
ists a characteristic radial position r x where the sign of the radial 
velocity v r changes: 

- Inward of r x , v r is negative, i.e., the flow is directed inward. 



1 The temperature at zo is sufficiently low that the matter is in molec- 
ular form at any radial position at the top of the atmosphere, and hence 
fi is close to 7/3 for the standard cosmic element mixture. 



- Outward of r x , v r is positive, i.e., the flow is directed out- 
ward. 

In model calculations of the solar nebula r x is typically located 
at a few or a few tens of AU and slowly changes with time (e.g. 
Ruden & Lin 1986, Paper II). This is far outside of the chemical 
active zone of the disk which is located inside of about 1 AU. 
Hence in one-zone models of the solar nebula, outward radial 
mixing of tracers from the chemical active zone to the outer 
solar system is only possible by turbulent diffusion (Paper II; 
Paper V). 

5.2. Meridional velocity field 

The real velocity field in disks is more complex than those of the 
one-zone approximation. Urpin ( 1984) considered higher order 
terms in the flow equations of the disk and found a characteristic 
flow pattern: 

- close to the disk midplane, the flow is directed outward 
whereas 

- at elevated heights above the midplane the flow is directed 
inward. 

This type of meridional flow field has been confirmed to ex- 
ist by 2-dimensional hydrodynamic calculations by Kley & 
Lin ( fl992l l, Rozyczk a et al . ( ff9"94l , Tscharnuter & Gail d2CT07l> . 
Regev & Gitelman (2002) conclude that the meridional flow 
field is a dynamical phenomenon that is of universal occurrence 
in accretion disksQ 

In the case of disks with a-viscosity, Takeuchi & Lin (2002) 
derived analytic expressions for the radial and azimuthal veloc- 
ity v r and v^, respectively. In their work v, was calculated in first 
order of the small perturbation parameter e ~ — ~ — . v^ is cal- 
culated in second order of e (which is the lowest non-vanishing 
order). Keller (2003) and Paper VI obtained the same result as 
Takeuchi & Lin (2002) and additionally obtained the vertical ve- 
locity v z in 1st order of s. 

In the following the meridional velocity field for disks with 
/3-viscosity is derived in the lowest non-vanishing order of s. 
According to Takeuchi & Lin ([2002b and Keller ( 120031 the ve- 
locity field in the stationary limit is given by 



v Ll , 



1 



m 



dP 



2GM t p dr 



vk 



d 



1 d , d v^ r d 

~T r P v 7> + ~TP V T V ^ 

rp or or r p oz oz 



2 

V'K 



i r d 

rp Jo or 



(41) 



(42) 



(43) 



As mentioned above, the lowest order, no-vanishing deviation of 
v,/, from Keplerian rotation is of second order in s. v r as well as 
v z are of first order in s. To simplify Eqs. ( f4Tb - j43l , the den- 
sity distribution of the stationary 1 + 1 -dimensional disk model is 
used, i.e. (e.g. Pringle [198TlPI 



p(r,z) - p c (r)e 



-r/(2H 2 ) 



(44) 



2 In the 3-dimensional extension the flow field of accretion disks 
should additionally show large vortices induced by the baroclinic in- 
stability (Klahr & Bodenheimer 2003 ; Klahr l2004l 

3 In Eq. ( 144b the self-gravity of the disk is neglected. It would be an 
interesting task to include the disk's self-gravity in the calculation of 
the meridional velocity field. 
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where 



Pc(r) 



\M\ ( R. 

6ttvH V r 



is the midplane density and 



H= J - 

2 Ok 



(45) 



(46) 



the scale height. R* denotes the radius of the protostar. Finally a 
radial power law for the midplane temperature is assumed, 



T c (r) oc r 



(47) 



where A is the radial temperature exponent! F° r these density 
and temperature profiles, the velocity field in yS-disks is as fol- 
lows: 



vAr,z) 



1 <r 1 7 l 

1_ (4 + A-/(r))-!-- -/- 

4 vi 4 r L 



I'K 



v r (r,z)=J3 


[>-l* 


v z (r,z)=/3 


H) 


with the factor 


™-(>/H 


and coefficients 


1 r 





i 1 

r/W - t(9 



O' 



vk 



V'K 



(48) 



(49) 



(50) 



d = - [(3 - i)(6 - 5A) - (18 - 8A)/(r)] 



C 2 = -(3-/0(9-5/1). 

To get a first impression of the meridional flow field see Fig. [4p. 
As in case of disks with a-viscosity (Takeuchi & Lin 2002, 
Paper VI) for the meridional velocity field of /3-disks (Eqs. (|48T > 
- (|5Q| >) there exists a characteristic height za where the radial 
velocity changes its sign: 



z ± _ 6-5A- 3/(r) 
H V 9-5/1 



(51) 



provided the quantity under the square root is positive, which 
either requires A < | (if fir) <K 1) or A > |. The latter is 
not realised in accretion disks models calculated with reasonable 
dust opacities. 

If A < I, which is satisfied in most parts of the disk (see 
Sect. 17.3. T| and Fig. |4jb)), v r is negative above za, i.e., the flow 
is directed inward, whereas v r is positive below za, i-e-, the flow 
is directed outward. Otherwise, in the case | < A < I, v r is 
negative at any height, and the disk matter (the mixture of gas 
and small particles) moves inward at all heights. 



4 In principle, this assumption is not needed. The quantity A in all 
equations for the velocity can simply be interpreted as the local value 
of A = -dlnT/dlnr. The impact of the vertical temperature structure 
on the meridional flow field would be also of interest. However, this is 
not considered in the present work. 



The factor f(r) affects the flow pattern only in the very inner 
parts of the disk (f(r — > oo) = 0). As a consequence the loca- 
tion of za> and therefore the flow pattern of the disk, essentially 
depends on the radial temperature exponent A. 

In the model calculations of Paper VI, A has been chosen as 
a fixed parameter and set constant throughout the disk. In the 
present work A is calculated locally, i.e., at each grid point from 
the radial temperature profile by a 3-point interpolation. The in- 
terpolated value of A is used to calculate the 2-dimensional ve- 
locity field of the disk (v r ,v z ) given by Eqs. j49l , < [50b - 

A vector transformation is implemented in the code for trans- 
forming the 2-dimensional velocity field (v r ,v z ) into polar coor- 
dinates (v' r ,Vg) since polar coordinates are used to calculate the 
transport of tracers (cf. Eq. ([8]l). 

Note int Eq. (l48l i that the angular velocity v^ deviates from 
the Keplerian velocity vk- In particular it varies with height z. 
However, the value of v^ from Eq. d48l > is not required in the 
present 2-dimensional model. 



6. Numerical treatment 

6. 1 . Initial and boundary conditions 

The boundary conditions are chosen as follows: 

- At the outer boundary (r = 200 AU) we set S = (no-torque 
condition). 

- At the inner boundary (r; = 0.096 AU) we choose the quasi- 
stationary boundary condition which has been introduced in 
Paper V. The quasi-stationary boundary condition removes 
the problem of the unphysical behaviour of "Z in the inner 
zone of the disk which appears in the case of the no-torque 
condition 2 = and ensures a smooth and physically more 
realistic radial distribution of the surface density (as well as 
other quantities) in the inner disk zone. 

- The concentrations cy at the outer boundary are chosen such 
that the silicates are amorphous (/ cr y,for = /cry.ens = 0) and all 
carbon that is not bound in CO is condensed into solid carbon 

(/car = 0.6). 

- The Ctj at the inner boundary are chosen such that the sil- 
icates are completely crystalline (/ C1 -y,f r = /cry.ens = 1) and 
solid carbon is completely oxidised (/ car = 0). 

- At the upper boundary of the 2-dimensional polar grid (#i = 
4.01°) the dust is assumed to have interstellar properties as is 
assumed for the outer boundary. Hence, the cy at the upper 
boundary are chosen such that the silicates are amorphous 
(/cry.for = /cry.ens = 0) and a fraction of / C ar = 0.6 of the total 
carbon is assumed to be bound in solid carbon grains, the 
remaining fraction being in CO. 

- For the upper boundary as well as the inner and outer bound- 
ary homogeneous von Neumann boundary conditions are 

dc- 

chosen for the tracres, i.e., we set -£■ — at the upper 

boundary and -S- - at both lateral boundaries. This choice 
prevents any flux of matter across the boundaries. 

- At the midplane of the disk symmetry conditions are applied. 



The initial conditions for the one-zone model are chosen as 
in Paper V for the model with solar abundance, i.e., the initial 
radial distributions of 2 and c,; car are calculated from a stationary 
model (dl,/dt = and dci iCm -/dt = 0) in which the silicates are 
assumed to be unaltered amorphous ISM grains. 

The cy at the position (r',6) of the 2-dimensional polar grid 
initially are set equal to the midplane value Cij(r = r'), i.e., the 
Cy initially are set constant in f?-direction. 
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Parameters (fi, M disk , y dlsk , . . .) 



T 



Initial model (S(r, t = 0), c u (r', 6, t = 0)) 



I 



Set of radial disk Eqs. im~i27l 


«« 




«« — i 






* 






NR iteration of (T c ,fi) up to A = 9- lfr 6 




1 




Solution of Eq. Il 12b for E(r) 




1 




Solution of the set of Eqs. d28t for the c,j(r) 




* 




New fayfttt(r), / cry , e „ s (/) and / car (r) 




1 




Fixpoint iteration of S(r) up to A = 10~ 5 








\ 




Set of vertical disk Eqs. d32b - l!36t 






«« — 1 






* 




Iterating up to an accuracy of A = 5 • 1CT 5 for za 




1 




Interpolation from the (r,z) grid onto the (r',9) grid 




1 




Solution of the set of Eqs. l[8j for the c, -,/(r', 0) 




1 




Interpolation from the (r',#) grid onto the (r,z) grid 




1 




New /c-y.forO, z), foty,em(r, z) and / car (r, z) 




1 




Fixpoint iteration of T(r', 6) up to A = 5 • 10" 3 








1 




New stellar mass M. from Eqs. d29t and d30b 




1 




Next time step 











Fig. 1. Flow chart of the code. The arrows to the right of the di- 
agram denote the iteration loops and the time loop, respectively. 
The inner two loops in the upper half refer to the one-zone model 
and the inner two loops in the lower half to the 2-dimensional 
model. A denotes the accuracy of the iteration. 'NR' is an ab- 
breviation for Newton-Raphson. The set of radial disk Eqs. ( [T3l 
- ( f27b is solved by a coupled Newton-Raphson method for the 
midplane temperature T c and the mean molecular weight /i. For 
details see text. 



6.2. 2-dimensional grids 

In the present model the disk structure is calculated in cylindri- 
cal coordinates (r,z) (1 + 1-dimensional approximation) whereas 
the tracer transport is calculated in polar coordinates (r',8). For 
a self-consistent computation of the tracer transport and the disk 
structure an interpolation of the relevant physical quantities be- 



tween both grids is applied. This is done by a 2D bilinear inter- 
polation (Press et al. 1992). 

The radial grid extends from 0.096 AU to 200 AU and con- 
sists of 332 logarithmically equidistant grid points (100 grid 
points per decade). The vertical grid size depends on the radial 
position and is variable in time since we adopt a step size con- 
trol for z in the calculation of the vertical disk structure. In the 
chemically active zone of the disk the number of vertical grid 
points usually is larger than 100 (maximum ~ 170 grid points) 
whereas it drops below 100 in the cool outer disk regions. The 
total number of grid points of the cylindrically symmetric grid 
adds up to about 35 000. 

For the polar grid an apex angle of G\ - 4.01° has been cho- 
sen. This choice ensures that the polar grid is embedded in the 
cylindrically symmetric grid at any time of the model calcula- 
tions. The polar angle G is discretised in L - 51 equidistant 
grid points 9/ between the upper boundary G\ - 4.01° and the 
midplane Gl - 0°. Note that the polar grid defines the 'domain 
of transport', i.e., the domain where the transport of tracers by 
advection and diffusion occurs during the model calculations. 
Beyond the polar grid no transport takes place (v r - v z - 0, 
D - 0). We did test calculations with L - 101 polar grid points 
as well as with an apex angle of G\ - 2.01° and found no sig- 
nificant deviations from the standard model with L — 51 and 
01 =4.01°. 



6.3. Numerical Solution 

The flow chart of the code of the present disk model is shown in 
Fig- HI The model calculations are performed as follows: 

First the radial disk structure in the one-zone approximation 
is calculated with standard methods (inner two loops in the upper 
half of Fig. [TJ see Sect.[3]and cf. Papers II and V). The equation 
for E, Eq. (fT2l) . is solved fully implicit. 

Secondly, the vertical disk structure is computed (innermost 
loop in the lower half of Fig. [TJ see Sect. |4|. In this way we 
obtain the disk structure in cylindrical coordinates (r,z). 

Thirdly, the 2-dimensional transport of tracers is calculated. 
This is done by 

1. Interpolating relevant quantities (T, D, n, «oh, v r , v z ) from 
the cylindrically symmetric grid (r,z) onto the polar grid 
(r',G), 

2. solving the set of 2D tracer equations (JHJ in polar coordinates 
(see Sect.[2jl, 

3. calculating the degrees of crystallisation of the silicates, 
/cry.for and / C ry,ens, as well as the degree of condensation of C 
in solid carbon, / cal , from the Cij(r', G) (see Eqs. (0 - (fTTTi). 
and 

4. interpolating / cr y,foi-, /oy.ens and / car from the polar grid (r', G) 
onto the cylindrically symmetric grid (r,z). 

To calculate the tracer transport self-consistently with the disk 
structure the temperature is iterated globally by a fixpoint itera- 
tion up to an accuracy of A g i = 5 ■ 10~ 3 (outer loop in the lower 
half of Fig. [TJ. The relative low accuracy of the global iteration 
is attributed to the errors which are introduced by the interpola- 
tion procedure. We did test calculations without performing the 
global iteration and found only small deviations from the stan- 
dard model with A g i = 5 ■ 10~ 3 which are not relevant for the 
final results. 

Finally, the stellar mass is updated (see Eqs. (|29l and ( TSOl l) 
before the next time step is executed (outermost loop in Fig.[TJ. 
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Table 1. Parameters used for the calculation of the disk models. 



Initial stellar mass 


M, fl 


\Mq 


Stellar effective temperature 


T, 


4 250 K 


Stellar luminosity 


L 


5L Q 


Stellar radius 


R, 


4.13 R Q 


Inner disk radius 


n 


5R, 
=0.096 AU 


Outer disk radius 


r a 


200 AU 


Molecular cloud temperature 


* cloud 


20 K 


Initial disk mass 


-Mdisk.O 


0.2 M Q 


Disk angular momentum 


-/disk 


10 53 gem 2 s -1 


Apex angle of the polar grid 


01 


4.01 ° 


Viscosity parameter 


P 


io- 5 


Number of polar 




332x51 


grid points 




= 16932 


Number of cylindrical 




-35 000 


grid points 







6.4. Model parameters 

The model parameters for the disk models calculated in the 
present work are shown in Table Q] The stellar and disk param- 
eters are chosen as in Paper V and are typical for solar like T 
Tauri stars and their surrounding accretion disks (see references 
given in Paper II). With respect to the choice of the value of the 
viscosity parameter /3 the reader also is referred to Paper V and 
the discussion therein. 

The model calculations are initiated with a small time step of 
10~ 8 yr which is increased slowly by an implemented time step 
control. The time step is limited to 50 yr to ensure numerical 
stability. We quit the model calculations at 10 6 yr. The model 
with the above standard parameters requires a CPU time of about 
two month on a P4 XEON 2.8 GHz computer. 



=2 

4 



0.01 



0.001 






0.01 



0.001 




r[AU] 

Fig. 2. Tracer transport in one-zone models with different ge- 
ometries at times t — (0), 10 5 (1) and 10 6 yr (2). The 'cylindri- 
cal' model IDC is shown with dashed lines and the 'polar' model 
1DP with solid lines, (a) Degree of crystallisation of forsterite 
/cry.for- (b) Fraction of C condensed in solid carbon / car . 



7. Results 

7.1. Models 1 DC and 1 DP 

Before we present the results of the 2D model calculations we 
first compare two model calculations in the one-zone approxi- 
mation which correspond to two different underlying disk ge- 
ometries. 

In the first model, henceforth called model IDC, the tracer 
transport is calculated in cylindrical coordinates (r,z) by assum- 
ing vanishing concentration gradients in the vertical direction, 

dc- 

i.e. , -rr = 0. The transport-reaction equation in the one-zone 
approximation then is given by Eq. d28l ). i.e., 



~ +Vr ~dr~ 



1 d n ddj Rij 
— — rnD —± + — - 
rn or or n 



(52) 



Model IDC essentially equals the standard model with solar el- 
ement mixture in Paper V. 

In contrast to model IDC, the tracer transport in the second 
model, henceforth called model 1DP, is calculated in polar co- 
ordinates {r',8) within a one-zone model. For this purpose we 
assume the concentration gradients in polar direction to be neg- 
ligible, i.e., -gjj- - 0. This assumption is almost identical to the 

assumption -|^ = made for Eq. (|52l . It is justified by the fact 
that the disk model is based on the idea of a flat and flaring disk 
(z <s r) which has a small apex angle. As well it is justified 
by the results of the 2D model calculations of the present work 



which show only slight concentration gradients in the vertical 
(respectively polar) direction since the vertical (polar) diffusion 
erases almost any vertical concentration gradient. 

With these basic settings the transport-reaction equation in 
2D polar coordinates, Eq. (0, transforms into 



dcu , dcu 

ot or 



1 d ,2 n dci 'J ^ Ri -J 
-^— — r nD — 1 

r' l n or' or' n 



(53) 



Equation d53l determines the tracer transport in the one-zone 
model 1DP 

Both Eqs. d52| > and d53l determine the tracer concentrations 
in the midplane of the disk where r — r' but the geometry of 
the disk in both models is different. In model IDC the disk is a 
flat slab whereas in model 1DP the disk resembles an outward 
flaring slab. For this reason in both models the diffusion terms in 
the transport-reaction equations, Eqs. d52l and ( f53l . differ from 
each other. 

The radial velocity profile in both models is assumed to be 
that of the one-zone model which is given by Eq. ((27}. 

Note that the radial density profiles n(r) that enters in the dif- 
fusion terms of Eqs. ( |53l and ( 1531 in both models is identically 
to the density profile of a flaring disk (cf. Eq. (fT7]i). However, 
to compare the tracer transport in models IDC and 1DP one ac- 
tually has to apply a density profile for model IDC with con- 
stant scale height to account for the cylindrical disk geometry in 
this model. In contrast, the polar disk geometry in model 1DP 
represents to be a good approximation of the real flaring disk 
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geometry. We realised this inconsistency after the termination 
of the model calculations. However, we refrain from recalculat- 
ing model IDC with a density profile with constant scale height 
since we intend to compare the model calculations of the present 
paper with those of the previous papers of this series which com- 
putes the tracer transport as model IDC, i.e., with cylindrical 
disk geometry and 'flaring' n(r). 

The dependence of the results for the tracer transport on the 
different assumed disk geometries is shown in Fig. [2] 

In Fig. [2^ the degree of crystallisation of forsterite / C ry,for 
at t = (0), 10 5 (1) and 10 6 yr (2) is plotted for the models 
IDC (dashed lines) and 1DP (solid lines). It clearly can be seen 
that the 'polar' model 1DP leads to a substantially less efficient 
radial mixing than the 'cylindrical' model IDC. In the region 
around 10 AU after 10 6 yr of the disk evolution, / cry ,for in model 
1DP is about one order of magnitude lower than in model IDC. 
This can be explained by the disk geometry. In model 1DP the 
disk flares. As a result in model 1DP tracers are diluted to a 
greater extent as they are mixed outward as compared to model 
IDC where the geometry is similar to a tube with parallel walls. 
However, polar coordinates obviously resemble more closely the 
real disk geometry than cylindrical coordinates do. Thus, in for- 
mer one-zone models of the solar nebula that treat the tracer 
transport in cylindrical geometry (Paper II; Paper V; Bockelee- 
Morvan et al. 2002) radial mixing has been overestimated. As a 
consequence the hypothesis of radial mixing that explains many 
interesting properties of primordial solar system bodies hardly 
can be maintained for previous one-zone models. For example 
the observed fraction of crystalline silicates of more than 10 % 
of the bulk silicate in many comets (e.g. Hanner et al. 1994, 
Crovisier et al. 1 19971 Hanner et al. 1997; Yanamandra-Fisher & 
Hanner [1999l Wooden et al . 120051120071 hardly can be explained 
by one-zone models with a realistic polar disk geometry. 

A way out of this dilemma is provided by 2-dimensional 
models with a realistic description of the flow field in the disk 
(see Sect. 1731). 

Carbon dust in model 1DP is mixed radially outward more 
efficiently than in model IDC (Fig.^). This can be explained 
by the dilution of the products of carbon combustion in the outer 
parts of the disk (e.g. CO, CH4; Gail 2002) which is more pro- 
nounced in the polar model 1DP than in the cylindrical model 
IDC. This again is due to the geometric effect described above 
for the radial mixing of the crystalline silicate grains. 



structures of enlarged grid point density are formed which stand 
out from the ambient grid with low resolution. An example is the 
arc-like structure that culminates at height z ~ 0.035 AU and ex- 
tends out tors 0.7 AU. This results from the transformation of 
enstatite to forsterite. Above and to the right of the arc enstatite 
and forsterite are both stable whereas beneath the arc enstatite is 
unstable and forsterite is the only stable silicate. The enstatite-to- 
forsterite transformation occurs at about 1 280 K. The two arcs 
culminating at z » 0.02 and 0.01 AU display the vapourisation 
fronts of forsterite (T ~ 1 360 K) and solid iron (T ~ 1 380 K), 
respectively. The vapourisation of corundum which is the most 
refractive dust species in the present model calculations occurs 
at r » 0.12 AU and T ~ 1 755 K (midplane values at t = 0). 

For comparison, in Fig. [3] also the radial run of the pressure 
scale height /z s (Eq. ( fT6l ); dashed line), the upper boundary of the 
2D polar grid (solid line) and the location where the atmosphere 
gets vertically optically thin (r = 1; big dots) is shown. The 
2D polar grid for tracer transport (not shown) is well embedded 
in the 1 + 1-dimensional cylindrically symmetric grid. Moreover, 
the 2D polar grid is well located within the vertically optically 
thick region of the disk. 

Note again that the 2D polar grid defines the area of turbu- 
lent and advective transport. Beyond these area no diffusive or 
advective transport is considered in the present model calcula- 
tions. 



7.3. 2-dimensional model calculations 

In the final step we present the results of the 2-dimensional 
model calculations. For the purpose of investigating the influ- 
ence of the disk's flow field on the tracer transport two model 
calculations with different flow patterns are performed: 

- In the first model, henceforth called model IDE, the flow 
field of the one-zone model (Eq. (1271 ) is extended to the disk 
regions off the midplane, i.e., at each polar grid point the po- 
lar radial velocity v' r is chosen equal to the one-zone radial 
velocity v r , v' r (r',ff) = v r (r - r'), and the polar velocity vg 
everywhere is set to zero, vg(r',6) = 0. Although this ve- 
locity field is 2-dimensional we call it in the following the 
'extended one-zone' velocity field (IDE). 

- In the second model, henceforth called model 2DM, the 
meridional velocity field is used which is calculated as de- 
scribed in Sect. [ 



7.2. 1+1 -dimensional grid 

In the next step we turn from the one-zone disk model to the 
1 + 1-dimensional disk model with 2-dimensional transport of 
tracers. 

To give an impression of the 1 + 1-dimensional cylindrically 
grid it is plotted in Fig. [3] at time t = of the standard model 
2DM for the innermost 0.8 AU of the disk and z > 0. As been 
described in Sect. 14.21 the solver for the disk vertical structure is 
equipped with an automatic step size control keeping the number 
of vertical grid points small. In particular the step size Az is small 
in such regions where physical quantities show large gradients. 
Such regions clearly can be seen in Fig. [3] 

On the one hand the grid point density gets high close to the 
upper boundary of the vertical grid as a consequence of the large 
pressure and density gradients. 

On the other hand the grid point density is increased in the 
narrow zones of dust vapourisation/destruction since tempera- 
ture and opacity are strongly variable there. As a result, arc-like 



7.3.1 . Flow fields in models 1 DE and 2DM 

To compare the extended one-zone velocity field and the merid- 
ional velocity field the flow fields of both models IDE and 2DM 
are plotted in Fig. [4] at time t — 10 5 yr in the region between 2 
andlOAU. 

The one-zone velocity field (Fig. 2^) is directed inward for 
most parts of the disk as it was mentioned in Sect. I5.ll The ve- 
locity increases inward as the viscous torque increases with de- 
creasing r. The absolute value of the velocity at 4AU, e.g., is 
20cms~ 1 whereas at the inner boundary it is 260cms~'. The 
feature with large grid point density at about 6 AU marks the ice 
condensation front. 

The meridional flow field (Fig. \4fa) shows a distinctly 
more complex structure than the one-zone velocity field. The 
meridional flow pattern shows large eddies which extend from 
both sides of the midplane of the disk to the disk's surfaces. 
Particularly the structure of the flow field in the region of the ice 
condensation front is of interest as the radial temperature distri- 
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Fig. 3. 1 + 1 -dimensional grid of the standard model at t — in the inner part of the disk. Each cross marks a grid point. The dashed 
line shows the radial profile of the pressure scale height h s (Eq. ( TTol l) and the solid line the upper boundary of the 2D polar grid. 
The big dots close to the upper boundary of the grid mark the location of the photosphere (t = 1). The arc-like structures mark the 
location of the vapourisation zones of strong absorbers (for details see text). 



bution shows a complex behaviour in that region. In detail the 
midplane radial velocity v' r (r', 8 = 0) shows the following radial 
dependence in model 2DM at t = 10 yr: 

(i) For r < 6 AU, v r is positive since for the radial temperature 

exponent in this region A < | holds (cf. Sect. l5.2t . 
(ii) In the range of 6 AU < r < 7.5 AU, v r is negative. There the 

opacity is a decreasing function of r that results in a radial 

temperature profile with A > | . 
(iii) In the range 7.5 AU < r < 9.5 AU, v r again is positive as the 

condensation front of ice leads to a temperature plateau, i.e., 

(iv) In the range 9.5 AU < 22 AU, v r for the same reasons as in 

(ii) is negative, 
(v) Finally, for r > 22 AU, v r is positive since the temperature 

adjusts to the ambient cloud temperature (r c i ou d = 20 K), i.e., 

A<l. 

Due to the accretion process the regions (i) - (v) slowly move 
radially inward. Additionally, the diffusive spreading of the disk 
causes a flattening of the radial temperature profile which nar- 
rows the ranges (ii) and (iv) with negative v r with time. Despite 
of this, the flow patterns (i) - (v) around the ice front persists for 
more than 10 6 yr of the disk evolution. The chemical active zone 
where in particular annealing of silicates and carbon combustion 
occurs, is always located within range (i). 

The polar velocity vg is more than one order of magnitude 
lower than v' r close to the midplane, and therefore plays no im- 
portant role with regard to mixing processes in the disk. 



7.3.2. Radial mixing: models 1DP and 1DE 

We now study the impact of the meridional flow field on the 
efficiency of radial mixing processes in the disk. For this pur- 
pose in Fig.[5]the midplane radial distributions (a) of the degree 
of crystallisation of forsterite / cry ,f i- and (b) of the fraction / car 
of C condensed into solid carbon are displayed for the models 
IDE (dashed lines) and 2DM (solid lines) with the same repre- 
sentation as in Fig. [2] For a comparison the result of the polar 
one-zone model 1DP is plotted in Fig. [5] with dotted lines. 

We first compare the models 1DP and IDE. In both mod- 
els the tracer transport is calculated in polar geometry with the 
same one-zone velocity field. Hence these models show the in- 
fluence on the tracer transport if extending it from ID to 2D. 
The main difference of models 1DP and IDE can be attributed 
to their different temperature structures. In the present work the 
models in the 1 + 1 -dimensional approximation generally lead to 
lower midplane temperatures for the chemically active parts of 
the disk than models in the one-zone approximation. This is an 
explicit result of the calculation of the vertical disk structure in 
the 1 + 1-models, particularly of calculating the radiative transfer 
in vertical direction. The midplane temperature difference be- 
tween the models 1DP and IDE can amount to more than 100 K 
in the innermost parts of the disk. As a consequence chemical 
reactions and dust processing in model IDE occur at a smaller 
distance r than in model 1DP, i.e., the annealing front of sili- 
cates (Fig. [5^) as well as the carbon combustion front (Fig.[5t>) 
in model IDE are shifted inward as compared to model 1DP. 
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Fig. 4. 2-dimensional flow field in the disk at t = 10 5 yr between 2 and 10 AU. The solid line marks the upper boundary of the polar 
grid (0i = 4.01°). (a) Model IDE. (b) Model 2DM. The arrows show the direction of the flow; their length is proportional to the 
velocity. For details see text. 



Superimposed to this temperature effect a pure 2D effect in- 
fluences the radial mixing. For evolution times t > 10 5 yr model 
IDE shows a larger concentration of crystalline silicates in the 
outer parts of the disk than model 1DP (Fig. [5^). This can be 
explained by the transition from 1- to 2-dimensional diffusion 
between both models. In model 1DP tracer particles experience 



pure radial random walks owing to the 1 -dimensional treatment 
of the diffusion. In contrast in model IDE tracers are radially 
mixed and additionally they are mixed in the vertical direc- 
tion, i.e., the tracer particles perform a 2-dimensional random 
walk. This leads to a more efficient radial mixing in model IDE 
than in model 1DP. In other words, concentration gradients are 
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Table 2. Degrees of crystallisation of forsterite and enstatite, f cly f or and f ciyfins , respectively, in the model calculations IDE and 
2DM at some selected times t and radial positions r. 





Model IDE 


Model 2DM 


/cry,for 


1AU 


5AU 


10 AU 


20 AU 


1AU 


5AU 


10 AU 


20 AU 


1.0- 10' yr 
6.5 ■ 10' yr 
1.0- 10 6 yr 


0.825 
0.219 
0.135 


0.056 
0.031 
0.021 


8.0 • 10~ 3 
0.017 
0.013 


1.9- 10- 4 
8.9 ■ 10~ 3 
8.6 • 10~ 3 


0.867 
0.624 
0.498 


0.248 
0.203 
0.134 


0.037 
0.075 
0.058 


4.6 ■ 10- 4 
0.016 
0.019 


/cry, ens 


1AU 


5AU 


10 AU 


20 AU 


1AU 


5AU 


10 AU 


20 AU 


1.0- 10' yr 
6.5 • 10' yr 
1.0- 10 6 yr 


0.658 
0.171 
0.103 


0.043 
0.024 
0.016 


5.9 • 10~ 3 
0.014 
0.010 


1.3- 10~ 4 
6.9 ■ 10~ 3 
6.7 ■ 10~ 3 


0.841 
0.599 

0.471 


0.231 
0.196 
0.128 


0.032 
0.073 
0.056 


3.6 ■ 10~ 4 
0.015 
0.018 



Table 3. Degree of condensation of carbon, f c[ 
positions r. 



in the model calculations IDE and 2DM at some selected times t and radial 





Model IDE 


Model 2DM 


/car 


0.5 AU 


1AU 


3AU 


5AU 


0.5 AU 


1AU 


3AU 


5AU 


1.0- 10' yr 
6.5 • 10' yr 
1.0- 10 s yr 


0.132 
0.499 
0.546 


0.441 
0.566 
0.582 


0.572 
0.591 
0.595 


0.589 
0.595 
0.598 


0.033 
0.188 
0.278 


0.132 

0.283 
0.368 


0.351 
0.430 
0.487 


0.473 
0.494 
0.571 



smoothed somewhat faster by 2-dimensional diffusion than by 
1 -dimensional diffusion. 



7.3.3. Radial mixing: models 1DE and2DM 

Finally we compare the extended one-zone flow field model IDE 
and the meridional flow field model 2DM. 

Figure[5]shows the remarkable impact of the meridional flow 
field on radial mixing. It clearly can be seen in the radial distri- 
bution of crystalline forsterite (Fig. [5^; enstatite yields a sim- 
ilar result). For example, after termination of the calculations 
at t - 10 6 yr the degrees of crystallisation of forsterite, / clT ,f r 
in model 2DM is up to 7 times larger than in model IDE (at 
rs3 AU). This result is a consequence of the meridional flow 
field for which the mass flux near the disk midplane is outward 
directed which makes radial mixing more efficient than in 1- 
dimensional models with pure influx. 

The region (ii) with negative v r at all heights (cf. Sect. 17.3.11 
and Fig. [4j}) only slightly impedes the outward mixing of crys- 
talline silicates in model 2DM because it is rather narrow. In this 
region outward transport is driven by diffusion only. In contrast, 
in region (iv) where v, < like in region (ii), the efficiency of 
outward mixing in model 2DM is significantly reduced because 
of the large extension of this region. At late stages of the disk 
evolution and for the very outer parts of the disk, the density of 
annealed silicates in model 2DM even falls below that of model 
IDE (at 10 6 yr for r > 35 AU). This is a consequence of the gen- 
eral outflow of matter in the one-zone velocity field of model 
IDE in the outer parts of the disk (cf. Sect. 15. It . The location of 
the transition from accretion flow to outflow in model IDE, r x , is 
located at 34 AU at time t — 10 6 yr. Note that r x moves outward 
during the disk evolution due to the disk spreading. 



In Sect. 17.11 we mentioned that the abundance of crystalline 
silicates in the comets hardly can be explained by model cal- 
culations that are based on a one-zone flow field and a polar 
disk geometry. Figure |6]reveals that this problem does not exist 
for the meridional flow field. The plot shows the time evolution 
of the degree of crystallisation of the silicates (forsterite: filled 
symbols; enstatite: open symbols) at the ice front, i.e., at the ra- 
dial position where water ice starts to freeze out onto grains (at 
T ~ 150 K). Model IDE is represented by squares and model 
2DM by circles. Note that the ice front slowly moves inward 
from 6.5 AU after t - 10 5 yr to 2.5 AU at t - 10 6 yr due to the 
accretion process. 

The degree of crystallisation of the silicates in model 2DM 
maintains values of above 10% beyond the ice front during 
the hole disk evolution. Note that the results prior to 10 5 yr 
are not physically meaningful since the disk has to evolve off 
from the initial model which is only a guess. The values of 
/cry.for and / C iy,ens exceed 20 % during most of the disk evolu- 
tion (t > 150.000 yr). In contrast, in model IDE the degree of 
crystallisation of the silicates never exceeds 10 % beyond the ice 
front where the comets must have been formed. Thus the abun- 
dance of crystalline silicates in many comets of more than 10 % 
can be explained by outward mixing of annealed grains from the 
inner disk parts into the region of comet formation. This is pos- 
sible only due to the existence of a meridional flow field. 

Note that crystalline silicates are absent in the interstellar 
medium (Kemper et al. 2004). Therefore it is improbable that 
the crystalline silicates in the comets originate from the ISM. 
This points to an origin of the crystalline silicates from the solar 
nebula. 

Figure \5]p displays the radial distribution of condensed car- 
bon in the same manner as Figure [5^. It can be seen that the 
narrow zone of carbon combustion in models IDE and 2DM 
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Fig. 5. Tracer transport in the 2-dimensional models IDE 
(dashed lines) and 2DM (solid lines) at times t - (0), 10 5 
(1) and 10 6 yr (2). Also shown is the result of the polar one-zone 
model 1DP (dashed dotted lines), (a) Degree of crystallisation 
of forsterite / cry? f or . (b) Fraction of C condensed in solid carbon 
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Fig. 6. Time evolution of the degrees of crystallisation of 
forsterite, / cry ,for, (filled symbols) and enstatite, /cry,ens, (open 
symbols) in the models IDE (squares) and 2DM (circles) at the 
position of the ice front. 



for different evolutionary periods is located at similar positions. 
This is a consequence of similar OH densities close to the car- 
bon combustion front in both type of models. However, owing 
to the meridional flow field, the carbon grains in model 2DM 
are diluted more efficiently than in model IDE. Hence in model 
2DM a substantial amount of the products of carbon combus- 
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Fig. 7. Cumulative plot of the dust-to-gas mass ratios /d ust : gas of 
the individual dust species versus radial distance in model 2DM 
after 5 x 10 5 yr of disk evolution. The individual dust species are 
labeled by obvious abbreviations. 



tion passes the ice front and is available for incorporation into 
the comets. This might explain in part the large abundances of 
hydrocarbons such as CH4 and C2H2 in comets (cf. Paper III). 

Tables |2] and [3] show values of / cry f or , / cr y,ens> an d /car for the 
models IDE and 2DM at some selected instants t and radial po- 
sitions r. They demonstrate quantitatively the results discussed 
above. 



7.3.4. Radial disk structure 

The radial disk structure is substantially modified by the merid- 
ional flow field. The opacity of the disk matter is reduced since 
carbon grains and amorphous silicate grains, that have a higher 
mass extinction coefficient than crystalline silicate grains, both 
are less abundant (see Fig. |5}. The lower opacity results in a 
more efficient cooling of the disk. As a result in the inner parts 
of the disk an extended radial zone establishes where the mid- 
plane temperature in model 2DM is more than 50 K below that 
of model IDE, e.g. from 0.4 to 2.3 AU. As the star accumulates 
matter from the disk this region moves slowly radially inward 
with time. As a result of the lower temperature the density is 
enhanced in model 2DM. As compared to model IDE the max- 
imum increase of the midplane density in model 2DM amounts 
to 15 %. 

The results show that it is essential for calculating the time 
dependent disk structure to consider both the most important 
opacity sources in the disk as well as their radial transport in 
the real flow field of the disk. 



7.3.5. Radial distribution of solids 

Figure [7] shows the radial distribution of solids in the disk for 
the meridional model 2DM at 5 x 10 5 yr as a cumulative plot 
of the dust-to-gas mass ratios of the individual dust species. 
From calculations of drag-induced drift of bodies in the so- 
lar nebula it is known that medium-sized bodies (diameter of 
~ 10 2 cm) located in the inner disk zone rapidly spiral into the 
protosun on timescales of about 10 4 yr (Weidenschilling et al. 
1989). Thus planetesimals have to be formed quickly from these 
bodies. Once formed, the planetesimals remain approximately 
in Keplerian orbits and mixing between different radii is not im- 
portant (Hayashi et al. 11985) . Hence plots of the kind of Fig. [7] 
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display the composition of the planetesimals at a given radial 
location if the time of their formation is known. 

Note that the headwind of the outward directed meridional 
flow field close to the midplane might lengthen the timescale of 
drag-induced lost in the sun and potentially speeds up the for- 
mation of planetesimals. 

8. Conclusions 

For a correct evaluation of the results of the model calculation 
we first consider the limitations of our present disk model. We 
list only the most important issues. 

First coagulation is not considered in the present disk model. 
As a consequence the dust particles remain small and are ver- 
tically mixed up to the surfaces of the disk by turbulent diffu- 
sion, and no sedimentation of grains toward the midplane takes 
place in the model. The disk structure (temperature, density etc.) 
will by strongly modified by coagulation owing to a change of 
opacity once this process commences. It would also be of great 
interest to know how big grains are radially transported by the 
meridional flow field. 

Secondly a more realistic treatment of the radiative transfer 
in the disk is important. This includes flaring by radiation from 
the protostar onto the disk surfaces, which is not considered in 
the model so far, as well as heat transport by convection. 

At present we apply the 1 + 1-dimensional approximation for 
calculating the disk structure while the tracer transport is com- 
puted in two dimensions. Hence, thirdly, exact 2D (or even 
3D) hydrodynamics would be desirable. Such calculations in 2D 
are currently underway and preliminary results are presented in 
Keller et al. d2004l > and Tscharnuter & Gail (12007b . 

Finally, in a complete model that extends the disk evolu- 
tion beyond the first million years and comprises the chemical 
and mineralogical evolution of growing bodies, the formation of 
gaps in the disk due to formation of massive planets should be in- 
cluded. It is to be expected that such gaps substantially influence 
or even suppress radial mixing of dust and gas species across the 
gaps. 

Although all these issues have to be adresses before we are 
able to predict, e.g., the mineral composition of the comets or 
the composition of planets, we draw the following preliminary 
conclusions: 

- The flow field of the disk substantially modifies the mix- 
ing of species in the disk. This is due to the fact that the 
timescales of mixing by advection is at least of same order 
as that for mixing by turbulent diffusion. Hence for tracking 
the spatial mixing of species, the real flow field of the disk 
has to be considered in model calculations of protoplanetary 
disks. In the present paper this is represented by the merid- 
ional flow field that already strongly enhances mixing com- 
pared to pure turbulent diffusion. It is to be expected that the 
effect will be even more enhanced in 3D-models. 

- The disk structure is modified by mixing processes. Species 
are chemically and morphologically alterated during the disk 
evolution, change their opacity properties and are mixed 
across the disk by diffusion and advection. This in turn af- 
fects the radiative cooling of the disk and, hence, the disk 
structure. To account for this, at least carbon combustion and 
silicate annealing have to be incorporated in realistic models 
of protoplanetary disks. The transport of these species has to 
be calculated self-consistently with the disk structure. 

- The tracer transport essentially depends on the geometry of 
the disk. Real disks flare, thus polar coordinates are more 



adequate for reproducing the disk geometry as compared to 
cylindrical coordinates. The flaring of the disk causes species 
to be significantly depleted in the outer parts of the disk as 
compared to the case of a flat geometry. Note that cylindrical 
coordinates have mostly been used in previous disk models 
that include tracer transport. 
- The depletion of tracers in models with polar disk geometry 
in the outer disk regions, particularly for annealed silicates, 
is more than compensated by the efficient outward mixing 
of these species via the meridional flow field. Therefore the 
presence of a fraction of crystalline silicates in some comets 
of more than 10% can easily be explained by the present 
model calculations. In contrast, the concentration of crys- 
talline silicates in the region of comet formation is more 
than an order of magnitude lower in models based on the 
one-zone flow field. This confirms also the results of the 2- 
dimensional model calculations of Keller & Gail (2004). The 
large abundance of methane and acetylene in the comets also 
could, at least in part, be due to the outward transport of ex- 
haustion gases by the meridional flow field. 

The results of the present model calculations have to be con- 
firmed by future hydrodynamic disk models. Fully 2-D radiation 
hydrodynamic models are presently in preparation, some pre- 
liminary results of which are already published in Tscharnuter 
& Gail ( 120071 1, but such calculations presently cannot compete 
with the long time basis and spatial resolution of model calcula- 
tions that can be obtained with 1 + 1 -dim models like the present 
one. 

As a next step, this model will be coupled with a detailed 
modeling of the chemistry in the disk. 
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